# M Orsi, 7 oct 09
# script to unwrap com xy coordinates of permeant from z-constraint simulation
# usage:  python /home/orsi/brahmsProject/analysisCode/latDiff/calcLatDiff.py

# parameters:
mTime = 0.1  # measurement time [mus] corresponds to 100 ns
stepDiff = 50 # this is the value set for the drug+hormone simulations [2009], might change!!!
timeStep = 20 * 10**(-9) # [microsecond = mus]
nValDiffuse = (int) (mTime/(stepDiff*timeStep))
#print nValDiffuse # should be 100000

inFile = open( 'xySoluteUnwrapped.dat', 'r' )

x = []; y = []   # lists 

for line in inFile:            # scan through file lines
    xval, yval = line.split()  
    x.append(float(xval)); y.append(float(yval)) # store values in x and y arrays

drx = x[len(x)-1] - x[0]
dry = y[len(y)-1] - y[0]
drSq = drx**2 + dry**2 # [nm^2]
print "Net displacement: %3.1f nm" % drSq**(0.5)

totSimulationTime = stepDiff * len(x) * timeStep  # [mus]
#print "Total simulation time: %g mus" % t
D = drSq / ( 4*totSimulationTime )
print "Lateral diffusion coefficient: %3.2f nm^2/mus" % D
outFile = open( 'xyDiffCoeff-200ns', 'w' )
outFile.write("%12.7f\n" % D);
drSq = drSq * (10**(-7))**2 # [cm^2]
#print "Net displacement: %3.2g cm" % drSq**(0.5)
totSimulationTime = totSimulationTime * 10**(-6) # [s]
#print "Total simulation time: %g s" % t
D = drSq / ( 4*totSimulationTime )
print "Lateral diffusion coefficient: %3.2e cm^2/s" % D

latDiff = 0;
for i in range(len(x)-nValDiffuse):
    drx = x[i+nValDiffuse] - x[i]
    dry = y[i+nValDiffuse] - y[i]
    drSq = drx**2 + dry**2 
    latDiff = latDiff + drSq / ( 4*mTime )

latDiff = latDiff / (len(x)-nValDiffuse)
print "Average lateral diffusion coefficient: %3.2f nm^2/mus" % latDiff
print "Average lateral diffusion coefficient: %3.2e cm^2/s" % (latDiff*10**-8)

outFile = open( 'xyDiffCoeff', 'w' )
outFile.write("%12.7f\n" % latDiff);

